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ABSTRACT 

In this paper we have examined numerically exact configurations of close binary sys¬ 
tems composed of incompressible fluids with internal flows. Component stars of bi¬ 
nary systems are assumed to be circularly orbiting each other but rotating nonsyn- 
chronously with the orbital motion, i.e. stars in binary systems have steady motions 
seen from a rotational frame of reference. We have computed several equilibrium se¬ 
quences by taking fully into account the tidal effect of Newtonian gravity without 
approximation. We consider two binary systems consisting of either 1) a point mass 
and a fluid star or 2) a fluid star and a fluid star. Each of them corresponds to gener¬ 
alization of the Roche-Riemann or the Darwin-Riemann problem, respectively. Our 
code can be applied to various types of incompressible binary systems with various 
mass ratios and spin as long as the vorticity is constant. 

We compare these equilibrium sequences of binaries with approximate solutions 
which were studied extensively by Lai, Rasio and Shapiro (LRS) as models for neu¬ 
tron star-neutron star (NS-NS) binary systems or black hole-neutron star (BH-NS) 
binary systems. Our results coincide qualitatively with those of LRS but are different 
from theirs for configurations with small separations. For these binary systems, our se¬ 
quences show that dynamical or secular instability sets in as the separation decreases. 
The quantitative errors of the ellipsoidal approximation amount to 2 ~ 10% for con¬ 
figurations near the instability point. Compared to the results of LRS, the separation 
of the stars at the point where the instability sets in is larger and correspondingly the 
orbital frequency at the critical state is smaller for our models. Since these sequences 
can be considered as evolutionary models of binary NS systems where component stars 
are approaching due to the effect of gravitational wave emission, we can expect that 
the final fate of such binary systems will be dynamical coalescence due to dynamical 
instability. 

Key words: binaries:close - hydrodynamics - instabilities - stars:neutron 
stars:rotation 


1 INTRODUCTION 


Close binary systems composed of massive compact objects 
such as neutron stars or black holes are relevant systems 
in astrophysics, because they are possible sources of many 
varieties of astrophysical phenomena. A close binary system 
of compact objects will evolve quasi-statically by decreasing 
its orbital separation due to gravitational radiation emis¬ 
sion until a final stage such as a violent merging or mass 
overflow will be reached. During inspiraling and coalescence 
phases, binary systems may be observed as sources of not 
only cosmological y-ray bursts (see e.g. Paczynski 1986) but 
also gravitational waves which will be detectable by kilome¬ 


ter size laser interferometers such as LIGO and VIRGO (see 
e.g. Abramovici et al. 1992; Thorne 1994). 


The evolution of bi nary systems prior to mergi ng wa s 
discussed by Kochanek (1992) and Bildsten & Cutler (1992). 
In their papers, it was pointed out that, for most viscosities 
known thus far, the time scale of angular momentum trans¬ 
fer within a neutron star is longer than that of evolution 
due to gravitational radiation reaction. Therefore the effect 
of viscosity may be negligible during the evolution of bi¬ 
nary systems toward coalescence. In other words, viscosity 
may not be powerful enough to synchronize spin motions 
nor change spin of component stars. 


In the framework of Newtonian gravity without viscos- 
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ity, circulation of a fluid remains constant even wh en the 
grav itat ional radiation reaction potential is included (Miller 
1974[^Jherefore, stars in binary systems will not rotate syn¬ 
chronously but have internal flows or spin in the frame of 
reference rotating with the orbital angular velocity. 


Because of a long time scale of energy dissipation due 
to gravitational radiation, such evolutionary sequences of 
binary systems can be followed by using equilibrium config- 
uratio ns even until stars come near to the innermos t stable 
orbit (see for example, Shapiro & Teukolsky 1983). How¬ 
ever, it is very difficult to construct stationary states of 
compressible stars with internal motions even today for non- 
axisymmetric configurations such as binary stars (Uryu & 
Eriguchi 1996, and references therein) except for the case of 
irrotational gaseous star binary systems (Uryu & Eriguchi 
1998a, 1998b). Several years ago, Lai Rasio & Shapiro (here¬ 
after LRS) developed a semi-analytical method to treat equi¬ 
librium of incompressible and compressible stars under the 
assumption that stellar shape s and cqui-density contours ar e 
approximated by ellipsoids ( Lai, Rasio fe Shapiro 1993a ) 
(LRS1 hereafter). In a series of papers, they applied their 
method to binary neutron star systems with spin and exten¬ 
sively examined the innermost stable circular orbits where 
dynamical instability of binary systems sets in (Lai, Rasio 
& Shapiro 1994a and 1994b: hereafter LRS2 and LRS3, re¬ 
spectively). 


Equilibrium configurations of incompressible binary 
fluid star systems with internal motions were computed nu¬ 
merically exactly by Eriguchi and Hachisu (1985). Further¬ 
more, Eriguchi, Futamase and Hachisu (1990) computed se¬ 
quences of binary disks with internal motions. In these pa¬ 
pers they computed sequences of equal mass binary disks or 
stars which can be used to study evolutionary sequences due 
to gravitational wave emission by keeping the circulation of 
each component constant and showed that there exist turn¬ 
ing points on equilibrium sequences. In other words, they 
indicated occurrences of secular and dynamical instabilities 
prior to the work of Lai, Rasio & Shapiro (LRS1). 


In this paper we have extended the computational 
method developed by Eriguchi and Hachisu (1985) so as 
to compute Roche type binary configurations and/or binary 
systems with arbitrary mass ratios and internal motions, 
and have reexamined equilibrium models of incompressible 
binaries whose parameters are in the physically meaningful 
range. First, we have computed sequences of synchronously 
rotating binary systems consisting of a point mass and a 
fluid, i.e. Roche type binary systems, and have compared our 
results with those obtained by using the ellipsoidal approxi¬ 
mation by LRS. Second, we have concentrated on configura¬ 
tions of irrotational Roche-Riemann and Darwin-Riemann 
problems for which vorticity seen from the inertial frame 
vanishes. We have compared our results again with ellip¬ 
soidal models computed in LRS1 and LRS2. Our results 
coincide qualitatively with those of LRS. For these binary 
systems, our sequences show that dynamical or secular in¬ 
stability sets in. The quantitative errors of the ellipsoidal 
approximation are usually 2 ~ 10% for models near the in¬ 
stability point. These sequences of equilibrium binary stars 
can be used as quasi-evolutionary models of binary NS sys¬ 
tems in which the separation of two component stars is de¬ 
creasing due to the effect of gravitational wave emission. 


Consequently such binary systems will result in coalescence 
due to dynamical instability. 

In the next section we present the basic equations and 
explain briefly the numerical method for incompressible fluid 
binary systems with spin, in section 3 we show computa¬ 
tional results of binary sequences for several types, and in 
section 4 we discuss evolution of binary systems of compact 
objects by applying our results. 


2 FORMULATION OF THE PROBLEM 
2.1 Assumptions and basic equations 


We assume plane symmetry of stellar structures about the 
x-y plane and the x-z plane in the rotational frame of refer¬ 
ence. Here the Cartesian coordinates ( x , y, z) are used. The 
2 -axis is chosen along the axis of rotation and centers of 
mass of two stars are located on the *-axis. In this paper we 
consider neutron stars which are approximated by incom¬ 
pressible fluids. Although this is a crude approximation, we 
have not found appropriate schemes to treat compressible 
self-gravitating gases with generic steady internal flows as 
mentioned in Introduction (see also Uryu & Eriguchi 1996). 
In particular, it is difficult to formulate a well-posed prob¬ 
lem and derive proper basic equations for steady states of 
compressible fluids, when the flow is parallel to the equato¬ 
rial plane (the planar flow) and the vorticity vector £ of the 
fluid (seen from the rotational frame of reference) is parallel 
to the rotational axis, £ = / e z where ( i = x, y, z) de¬ 
notes the unit basis vector along each coordinate axis. These 
assumptions are not only physically natural, since the grav¬ 
itational radiation reaction potential conserves the circula¬ 
tion of the flow field (see below), but also necessary to make 
the problem tractable, since it reduces the number of phys¬ 
ical variables. However, for compressible models, we cannot 
write down the basic equations which are consistent with 
the number of independent physical variables. This occurs 
because the flow becomes independent of 2 under the above 
assumptions Q (see also the formulation for incompressible 
case below). 

On the other hand, for incompressible fluids we can 
formulate the problem properly by introducing the stream 
function for planar flows as follows: 


U x 


<9T 

~dp 


Uy 


< 94 > 

dx 


(1) 


where and ( u x , u y ) are the stream function and x- and 
^-components of the velocity vector in the rotational frame, 
respectively. Note that 4/ and (u x , u y ) depend only on (x, y) 
and that the equation of continuity is automatically satis¬ 
fied. Thus we can write down basic equations and boundary 
conditions consistently. 

If the stream function is a function of the 2 -component 
of vorticity vector, £, which we will refer to as the vortic¬ 
ity, the Euler equation of fluid motion can be integrated 


* However, we can formulate and solve the problem of compress¬ 
ible fluid binary consistently for irrotational case, i.e. the vortic- 
i tv of fluid vanishes iden tically in the inertial frame of reference 
( Uryu fe Eriguchi 1998h| ). In this case, we do not need to require 
the planar flow explicitly. 


© 0000 RAS, MNRAS 000, 000-000 








Incompressible fluid binary 3 


to the Bernoulli equation. Since we consider the case with 
C, = const everywhere in the fluid, the Bernoulli equation 
becomes as follows: 

\{ul + ul) + (C + 2n)* 

-in 2 (X 2 +y 2 )+ - + </> = C, (2) 

2 P 

where fl, p, p, <j> and C are the orbital angular velocity, the 
pressure, the density, the Newtonian gravitational potential 
and the Bernoulli constant, respectively. The quantity X 2 + 
Y 2 means the square of the distance of fluid element from 
the rotational axis. 

By substituting the velocity components Eq. (1) into 
the following definition of vorticity 

C= dUy_dU* 

ox oy 

the equation for the stream function can be derived as, 

A® = -C, (4) 


where 


A = 


dx 2 dy 2 


( 5 ) 


As for the gravitational potential, we use the integral 
form of the Poisson equation as follows: 


4>(r)=-GpJ^ | r J 1 r , | d 3 r . (6) 

Here G is the gravitational constant and V is the total vol¬ 
ume over which the fluid extends. 

Boundary conditions for the equations described above 
are as follows: 


angle measured from the positive direction of the x-axis. Re¬ 
lations between the Cartesian coordinates (X±,Y±, Z±) and 
the polar coordinate systems (r±,9±, p±) are expressed as 
follows: 

X± = ±d± + r± sin 9± cos p± , 

Y± = ±d± + r± sin#± sini£± , (10) 

Z± = r± cos 9± . 

Hereafter we will omit subscripts ± except for necessary 
cases in order to make expressions simple. Since deformation 
of the stellar surface from sphere is not so large for binary 
stars, we may consider the stellar surface as a function of 9 
and p and write it as R(9, p) in this coordinate. 

In the coordinate systems adopted here, we can ana¬ 
lytically integrate the gravitational potential (sj) along the 
r-direction from the center to the surface of the star R(9, ip) 
because the density p is constant. For incompressible stars, 
we only need to solve for the shape of the surface R{9, p) 
because physical quantities such as the pressure or the gravi¬ 
tational potential are completely determined from the shape 
of the surface (see e.g. Eriguchi, Hachisu & Sugimoto 1982; 
Eriguchi & Hachisu 1983; Eriguchi & Hachisu 1985). 

The stream function T is also expressed in terms of the 
surface variable. Since £ is assumed to be constant, we can 
explicitly write down the stream function by using power 
series expansion on the equatorial plane as follows: 

flt(r sin 8, p) = 

{ OO -v 

i(r sin 9) 2 + E a m (r sin 9) m cos m p > . (11) 

m=o ' 

Imposing the boundary condition (^), i.e. 


p(x,y, z) = 0, on the surface of the star, 


(7) *(R(n/2,p),p) = 0, 


( 12 ) 


and 

V(x,y) = A, 

along the surface on the equator of the star, (8) 

where A is a constant. Since we can choose its value ar¬ 
bitrary, we set A = 0 for convenience. It should be noted 
that boundary conditions for the gravitational potential are 
already included in the integral representation Eq. (^). 

2.2 Solving method 

In actual computations, we introduce two spherical coor¬ 
dinate systems (r±,0±,<^±) whose origins are fixed on the 
x-axis at distances 

_ I ^°± ut + ^ I 

d± = - - - , (9) 

for each fluid component, where R™* and R± are distances 
from the rotational axis to the outer and inner edges of 
the star on the x-axis, respectively. Here subscripts + and 
— correspond to the primary and the secondary stars, re¬ 
spectively. For Roche type binary systems, the primary star 
corresponds to a fluid star and the secondary to a point 
mass. Therefore note that usually the mass of primary M+ is 
smaller than that of secondary M_, i.e. M+ < JWT. The an¬ 
gle 8± is the zenithal angle measured from the positive direc¬ 
tion parallel to the 2 -axis and the angle p± is the azimuthal 


on the equatorial surface, we can obtain coefficients a m for 
a given vorticity. Components of the fluid velocity in the ro¬ 
tational frame (u x ,u y ) are computed by differentiating this 
expression. 

Substituting the above formulae into the Bernoulli 
equation on the fluid surface and imposing the boundary 
condition (|t[) for the pressure there, we finally obtain the 
basic equations for numerical computations. They are basi¬ 
cally the equations (^) and (|l^) (see Appendix A). The de¬ 
rived equations contain two unknown variables R(9, tp) and 
am, and several unknown physical constants. These equa¬ 
tions are calculated in the coordinate system (9 , p) over the 
region 9 £ [0,7r/2] and p £ [0,7r]. This region is discretized 
equidistantly into grid points as follows: (0; e , pi v ) (0 < ig < 
Ng, 0 < i v < N v ) and we choose (Ng,N v ) = (16,32). As 
for the coefficients a m , we include terms up to m = N v / 2. 
The trapezoidal formula is used for the integration of the 
gravitational potential. 

One model of a binary system can be numerically com¬ 
puted by specifying some physical quantities which charac¬ 
terize the binary system. We can choose the values of phys¬ 
ically meaningful parameters so as to be able to compute 
an appropriate sequence by changing some of them. For the 
present purpose, it is convenient to specify the mass ratio 
and the separation of two component stars as two parame¬ 
ters. For other parameters, we choose relations for the vor¬ 
ticity of the fluid star (see Appendix B). Since we consider 


© 0000 RAS, MNRAS 000, 000-000 




4 K. Uryu and Y. Eriguchi 


equilibrium sequences whose vorticity is conserved in the 
inertial frame as stars evolve adiabatically by the force of 
the gravitational radiation reaction potential, the fol lowing 
qua ntity must be constant in time as well as in space (Miller 
1974F1 

£ + 2 12 = £o = constant, 


(13) 


where Co is the vorticity in the inertial frame of reference. 
In particular, the stars with Co = 0 are usually referred to 
as irrotational configurations. 

The basic equations are discretized on the mesh points. 
Since these difference equations are nonlinear equations for 
R(9, if), a m and unknown constants, C and 12, the Newton- 
Raphson method is applied to solve them. Since the Newton- 
Raphson scheme requires a rather large matrix, we were able 
to employ only the mesh number mentioned above. For typ¬ 
ical models, our computational time was ~ 3.5 minutes for 
1 iteration for a fluid star-fluid star binary system and ~ 5.7 
seconds for a point mass-fluid star binary system by using 
a Fujitsu VX vector computer. A further discussion on the 
solving method is described in the appendices in detail. 


2.3 Consistency at the stellar surface boundary 


The formulation of the problem mentioned above is mathe¬ 
matically well-defined as a boundary value problem for the 
elliptic partial differential equation. Except for rigidly rotat¬ 
ing models, this formulation is essentially the same as the 
classical problem of two-dimensional vortex flow in a rotat¬ 
ing frame. We implicitly solve the equation of continuity by 
introducing the function and hence we treat all equations 
consistently. Actually the problem can be solved numerically 
by using the above prescription whose results we will show 
below. 

However, at the present time, we cannot show analyti¬ 
cally that the stream function T on the equatorial plane is 
consistent with the boundary condition on the stellar sur¬ 
face. We compute the stream function on the equatorial 
plane. On the other hand, the stellar surface is determined 
by the condition of equation (^). Since the flow lines corre¬ 
spond to the 4/ = constant lines, these should coincide with 
the p = constant lines. In other words, since the velocity vec¬ 
tor is always parallel to the stream line, the velocity vector 
in the rotating frame u and the normal vector of the stellar 
surface n are orthogonal everywhere on the surface. In order 
to check our computational results, we have computed the 
following quantities on the stellar surface numerically: 


<i = 


(14) 


To be real solutions, this quantity should vanish on the stel¬ 
lar surface. The normal vector n can be written as follows 
in the spherical coordinates: 


1 OR _ 1 dR 

R(8, ip) 89 e ° R(9, p) sin 9 8 <p e<p 


(15) 


where e r , eg and e v are the the unit basis vectors of the 
spherical coordinates. We have computed the quantity q nu¬ 
merically at every grid point on the discretized stellar sur¬ 
face. Although we have carried out numerical differentiation 
to obtain this quantity, the value is always less than 0.05, i.e. 
q < 0.05 everywhere on the stellar surface even for highly 


deformed configurations. Therefore, our method gives con¬ 
sistent solutions for stellar configurations as we expect. One 
can also find an excellent coincidence of the present results 
with other results which have been computed by using a to¬ 
tally different formulation and a numerical method (Uryu & 
Eriguchi 1998a, 1998b). This also ensures the validity of our 
present results. 


3 RESULTS 


It is convenient to use non-dimensional quantities to com¬ 
pare our results with those of others. Thus physical quanti¬ 
ties are normalized by using normalization factors as follows. 

The non-dimensional orbital angular velocity Cl is de¬ 
fined as: 


12 = 


12 

YG~P 


(16) 


Concerning the separation, the non-dimensional separation 
between centers of mass of two component stars do is defined 
as 


( Ag+ + d,G -) 

(3M tot /4npy/ 3 ’ 


(17) 


where the total mass is expressed as Mtot = M+ + A/_, 
and we define d,G± as follows: 


Ag± 


p\J v X ± dV\ 

M± 


(18) 


The separation between two geometrical centers is normal¬ 
ized as 


(d+ + d—) 

(3 Mtot/^np) 1 / 3 


(19) 


The total angular momentum J and the total energy E are 
normalized as 


YgmU?p -VO ’ 


( 20 ) 


and 


E = 


GM^pG 3 ' 


( 21 ) 


In the Tables we also show the luminosity of gravita¬ 
tional radiation, which is related to the quadrupole moment 
tensors of the binary system. In the rotating frame of refer¬ 
ence, the quadrupole moment tensors are defined as 


r(rot) 

ij 


= P / (xiXj - 8i 


-)dV, 


( 22 ) 


where i, j = 1,2,3 correspond to the x, y, ^-component of the 
Cartesian coordinates, respectively. By taking the symmetry 
of the system into account, the quadrupole formula for the 
energy generation rate is written as follows: 

dE _ 32 G 0 6 t T (rot) r(rot)^ 2 roo\ 

dt ~ 5 c 5 1 11 22 ’ ’ [ ’ 

where c is the velocity of light. Derivation of this formula 
is shown in the book of Misner, Thorne & Wheeler (1970) 
(see also Eriguchi, Futamase & Hachisu 1990). This energy 
generation rate is normalized as 
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dE_ 

dt 



(24) 


Table 1. Equilibrium sequence of Roche type binary systems 
with equal mass 


where non-dimensional time is defined as 

t = Cfpt , (25) 

and 

vo = yfcM^y/ 3. (26) 

In the following subsections, our results for three differ¬ 
ent types of binary configurations are discussed. They are, 

(i) Roche configurations with M+/M- = 1 and 0.1 . 

(ii) Irrotational Roche-Riemann configurations with sev¬ 
eral mass ratios. 

(iii) Irrotational Darwin-Riemann configurations with 
several mass ratios. 

LRS computed several types of binary configurations by us¬ 
ing the ellipsoidal approximation. Concerning binary solu¬ 
tions, nothing but synchronously rotating fluid fluid binary 
configurations with equal mass, i.e. Darwin models, are com¬ 
pared in their papers with fully deformed solutions. We com¬ 
pare our results with those of LRS for models listed above 
and clarify quantitative differences. These comparisons pro¬ 
vide the evidence that our computational method works ac¬ 
curately. 


3.1 Roche type binary systems 


We show our results in Tables [1|, |2| and |3| for Roche sequences 
with the mass ratio M+/M- = 1,0.5 and 0.1, respectively. 
These sequences can be regarded as simplified mod els of BH 
NS b inary systems in the limit of strong viscosity (Kochanek 


1992 ) and in the frame of Newtonian gravity. 

In Figures [j] and we show our results and those of 
LRS1. In Figures jl](a)-(c) and !( a ) -(c), physical quanti¬ 
ties, J, E and Cl, are plotted against the separation of mass 
centers of two component stars dG and against the non- 
dimensional angular momentum. The models on a sequence 
at the smallest separation in Figures fia), (b) and l a ), (b) 
almost correspond to those at the Roche limit. Here we de¬ 
fine the Roche limit as the smallest separation of two centers 
of mass of two bodies. Note that models shown in the last 
entry of the Tables are not the exact critical models. 

As seen from these figures, two results obtained by two 
different methods agree rather well even near the turning 
points of sequences. Differences of J or I? between two re¬ 
sults are less than 0.5% for most part of the sequence. It 
is also shown that our values of the separation between the 
centers of mass of the fluid star and the point mass agree 
with those of LRS to within 2 ~ 6% around the models 
with the minimum of J or E where secular instability sets 
in and also around the models near the Roche limit stage. 
Our results show that the separations where the instability 
sets in increase by the amount of roughly ~ 4% and ~ 2% 
for M+/M- = 1 and 0.1 sequences, respectively, and cor¬ 
respondingly the orbital frequency at the critical state de¬ 
creases by the amount of 6% and 3%, respectively, compared 
to the results of LRS. 

Here it should be noted that the secular instability oc¬ 
curs via angular momentum transfer between the orbital one 


d 

da 

a 

J 

E 

dE/dt 

5.0 

4.079 

2.486(-l) 

4.104(-1) 

-3.525(-l) 

3.883(-3) 

4.8 

3.930 

2.630(-l) 

4.039(H) 

-3.542(-l) 

4.688(-3) 

4.6 

3.782 

2.787(-l) 

3.974(H) 

-3.560(-l) 

5.693(-3) 

4.4 

3.635 

2.959(-l) 

3.909(H) 

-3.579(-l) 

6.957(-3) 

4.2 

3.490 

3.146(-1) 

3.845(H) 

-3.598(-l) 

8.556(-3) 

4.0 

3.347 

3.351(-1) 

3.782(H) 

-3.619(-1) 

1.059(-2) 

3.8 

3.207 

3.576(-l) 

3.720(H) 

-3.640(-l) 

1.319(-2) 

3.6 

3.070 

3.821(-1) 

3.661(H) 

-3.662(-l) 

1.653(-2) 

3.4 

2.937 

4.088(-l) 

3.606(H) 

-3.684(-l) 

2.083(-2) 

3.2 

2.810 

4.377(-l) 

3.555(H) 

-3.706(-l) 

2.637(-2) 

3.0 

2.689 

4.686(-l) 

3.510(H) 

-3.726(-l) 

3.347(-2) 

2.9 

2.632 

4.847(-l) 

3.491(-1) 

-3.735(-l) 

3.773(-2) 

2.8 

2.577 

5.010(-1) 

3.475(-l) 

-3.744(-l) 

4.249(-2) 

2.7 

2.526 

5.175(H) 

3.461(-1) 

-3.751(-1) 

4.780(-2) 

2.6 

2.478 

5.339(-l) 

3.451(-1) 

-3.756(-l) 

5.366(-2) 

2.5 

2.434 

5.500(-l) 

3.446(-l) 

-3.759(-l) 

6.004(-2) 

2.4 

2.395 

5.653(-l) 

3.445(-l) 

-3.759(-l) 

6.685(-2) 

2.3 

2.362 

5.794(H) 

3.450(-l) 

-3.757(-l) 

7.387(-2) 

2.2 

2.337 

5.916(H) 

3.460(-l) 

-3.751(-1) 

8.076(-2) 

2.1 

2.320 

6.010(H) 

3.477(-l) 

-3.741(-1) 

8.690(-2) 

2.0 

2.312 

6.065(H) 

3.495(-l) 

-3.730(-l) 

9.131(-2) 

Table 2. Equilibrium t 

sequence of Roche type binary systems 

with 

M+/M- 

=0.5 




d 

da 

h 

3 

E 

dE/dt 

5.0 

3.643 

2.948(-l) 

3.422(-l) 

-2.031(-1) 

5.427(-3) 

4.8 

3.519 

3.107(H) 

3.370(-l) 

-2.047(-l) 

6.472(-3) 

4.6 

3.397 

3.277(H) 

3.318(-1) 

-2.063(-l) 

7.749(-3) 

4.4 

3.277 

3.461 (-1) 

3.268(-l) 

-2.080(-l) 

9.312(-3) 

4.2 

3.159 

3.658(H) 

3.218(-1) 

-2.098(-l) 

1.123(-2) 

4.0 

3.045 

3.869(H) 

3.170(-1) 

-2.116(-1) 

1.359(-2) 

3.8 

2.935 

4.093(H) 

3.125(-1) 

-2.134(-1) 

1.648(-2) 

3.6 

2.830 

4.330(H) 

3.082(-l) 

-2.152(-1) 

2.000(-2) 

3.4 

2.730 

4.577(H) 

3.044(-l) 

-2.169(-1) 

2.427(-2) 

3.2 

2.638 

4.830(H) 

3.011(-1) 

-2.185(-1) 

2.934(-2) 

3.1 

2.596 

4.956(H) 

2.997(-l) 

-2.192(-1) 

3.218(-2) 

3.0 

2.556 

5.081(H) 

2.986(-l) 

-2.198(-1) 

3.521 (-2) 

2.9 

2.520 

5.201(H) 

2.976(-l) 

-2.202(-l) 

3.840(-2) 

2.8 

2.487 

5.316(H) 

2.970(-l) 

-2.206(-l) 

4.170(-2) 

2.7 

2.459 

5.422(H) 

2.967(-l) 

-2.207(-l) 

4.501(-2) 

2.6 

2.435 

5.516(H) 

2.967(-l) 

-2.207(-l) 

4.823(-2) 

2.5 

2.418 

5.593(H) 

2.972(-l) 

-2.205(-l) 

5.115(-2) 

2.4 

2.407 

5.649(H) 

2.979(-l) 

-2.201(-1) 

5.352(-2) 

2.3 

2.403 

5.677(H) 

2.989(-l) 

-2.195(-1) 

5.501(-2) 

2.2 

2.405 

5.676(H) 

2.996(-l) 

-2.191(-1) 

5.527(-2) 


and the spin because of the existence of viscosity (see e.g. 
Aizenman 1968; Hachisu 1979; Hachisu & Eriguchi 1984a, 
1984b, 1985; Lai, Rasio & Shapiro 1993a). Consequently, as 
the separation is decreased, the system reaches first a state 
where secular instability sets in and next comes to the Roche 
limit state. This behavior is also seen in models of LRS1 for 
sequences whose mass ratios are the same as those of ours. 


3.2 Irrotational Roche—Riemann type binary 

systems 

As mentioned in Introduction, neutron stars may be com¬ 
posed of inviscid fluids (Kochanek 1992; Bildsten & Cutler 
1992). For incompressible fluids, the vorticity Co is conserved 
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Table 3. Equilibrium sequence of Roche type binary systems 
with M+/M—=0.1 


d 

da 

n 

J 

E 

dE/dt 

10.0 

4.630 

2.055(-l) 

1.407(-1) 

-3.215(-2) 

2.239(-4) 

9.6 

4.461 

2.173(-1) 

1.382(-1) 

-3.269(-2) 

2.699(-4) 

9.2 

4.293 

2.302(-l) 

1.356(-1) 

-3.326(-2) 

3.271(-4) 

8.8 

4.127 

2.442(-l) 

1.330(-1) 

-3.388(-2) 

3.988(-4) 

8.4 

3.963 

2.596(-l) 

1.304(-1) 

-3.454(-2) 

4.890(-4) 

8.0 

3.800 

2.764(-l) 

1.278(-1) 

-3.524(-2) 

6.032(-4) 

7.6 

3.641 

2.949(-l) 

1.252(-1) 

-3.599(-2) 

7.485(-4) 

7.2 

3.485 

3.150(-1) 

1.226(-1) 

-3.678(-2) 

9.336(-4) 

6.8 

3.332 

3.370(-l) 

1.200(-1) 

-3.762(-2) 

1.170 (-3) 

6.4 

3.185 

3.608(-l) 

1.175(-1) 

-3.851(-2) 

1.473(-3) 

6.0 

3.043 

3.865(-l) 

1.150(-1) 

-3.943(-2) 

1.857(-3) 

5.6 

2.909 

4.139(-1) 

1.127(-1) 

-4.037(-2) 

2.341(-3) 

5.2 

2.784 

4.425(-l) 

1.105(-1) 

-4.130(-2) 

2.937(-3) 

4.8 

2.673 

4.712(-1) 

1.086(-1) 

-4.217(-2) 

3.646(-3) 

4.4 

2.579 

4.984(-l) 

1.071(-1) 

-4.290(-2) 

4.433(-3) 

4.3 

2.559 

5.046(-l) 

1.068(-1) 

-4.305(-2) 

4.631(-3) 

4.2 

2.540 

5.104(-1) 

1.066(-1) 

-4.318(-2) 

4.827(-3) 

4.1 

2.524 

5.158(-1) 

1.064(-1) 

-4.328(-2) 

5.014(-3) 

4.0 

2.510 

5.207(-l) 

1.062(-1) 

-4.336(-2) 

5.192(-3) 

3.9 

2.497 

5.249(-l) 

1.061(-1) 

-4.342(-2) 

5.354(-3) 

3.8 

2.488 

5.284(-l) 

1.061(-1) 

-4.344(-2) 

5.495(-3) 

3.7 

2.481 

5.311(-1) 

1.061(-1) 

-4.343(-2) 

5.610(-3) 

3.6 

2.478 

5.328(-l) 

1.062(-1) 

-4.339(-2) 

5.693(-3) 

3.5 

2.477 

5.334(-l) 

1.063(-1) 

-4.331(-2) 

5.737(-3) 

3.4 

2.480 

5.330(-l) 

1.065(-1) 

-4.322(-2) 

5.740(-3) 

3.3 

2.484 

5.321(-1) 

1.066(-1) 

-4.314(-2) 

5.715(-3) 


during evolu tion due to the gravitational radiation reaction 
( |Miller 1974 ). Furthermore the orbital angular velocity of 
the star at a large separation is much smaller than that 
at a separation of the stellar size. Since £o is of the order 
of the orbital angular velocity at large distances, the con¬ 
stant in equation (|I^) can be regarded very small and can 
be approximated to zero. In other words, almost all real¬ 
istic neutron star binary systems can be approximated by 
irrotational configurations as far as viscosity is small. 

We show our results for irrotational Roche-Riemann se¬ 
quences with M+/M- = 1,0.5 and 0.1 in Tables SI and 
P, respectively. In Figures y and [i], our results and those 
of LRS1 are displayed for the sequences with M+/M- = 1 
and 0.1. For all the sequences, the turning points of J and E 
against the separation exist as seen in Figures and []|, which 
is qualitatively the same as that of the ellipsoidal approxi¬ 
mation. It implies that there is a point where the dynamical 
instability will set in on the sequence between large separa¬ 
tion states and the Roche-Riemann limit of the sequence. 
Here the Roche-Riemann limit is defined as the generaliza¬ 
tion of the Roche limit to fluid stars with internal flows. 
Therefore, it is probable that an irrotational star around 
a massive point source will fall onto the point mass due 
to dynamical instability as the separation becomes smaller. 
This behavior is the same as that obtained by the ellipsoidal 
approximation. Under the ellipsoidal approximation, there 
always exists a dynamical instability point for any mass ra¬ 
tio for Roche-Riemann type binary systems before reach¬ 
ing the Roche-Riemann limit as the separation decreases 
(LRS1). The quantitative difference of the physical values 
between these two results near the turning point of the se¬ 
quence amounts roughly to 3 ~ 5% as seen from Tables and 
Figures. In particular, the separation of the minimum of J 


Figure 1. Physical quantities of a sequence of Roche type bi¬ 
nary systems with equal mass, (a) Total angular momentum as a 
function of a binary separation, (b) Total energy as a function of 
a binary separation, (c) Orbital angular velocity as a function of 
the total angular momentum. Dashed and solid curves denote the 
sequence of LRS1 and our sequence, respectively. See text about 
the normalization factors for each quantity. 


and E increases by 3.5% and 2.5% for the sequences with 
M+/M- = 1 and 0.1, respectively, and the orbital frequency 
decreases by 5% and 4%, respectively. 


3.3 Irrotational Darwin—Riemann type binary 
systems 

In LRS2 they calculated physical quantities for irrotational 
Darwin-Riemann sequences by assuming ellipsoidal config¬ 
urations. As shown in the previous sections, exact compu¬ 
tations have given qualitatively the same results for models 
even with smaller separations. Thus we have further com¬ 
puted stationary solutions of two fluid stars with (jo = 0. 
These types of binary systems can be regarded as mod¬ 
els of viscosity-free NS-NS binary systems just before co- 
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Table 4. Equilibrium sequence of irrotational Roche—Riemann 
type binary systems with equal mass 


d 

dG 

n 

J 

E 

dE/dt 

5.0 

4.057 

2.506(-l) 

3.969(-l) 

-3.544(-l) 

3.987(-3) 

4.8 

3.907 

2.653(-l) 

3.895(-l) 

-3.563(-l) 

4.826(-3) 

4.6 

3.757 

2.814(-1) 

3.821(-1) 

-3.583(-l) 

5.877(-3) 

4.4 

3.609 

2.990(-l) 

3.746(-l) 

-3.605(-l) 

7.204(-3) 

4.2 

3.463 

3.182(-1) 

3.670(-l) 

-3.629(-l) 

8.886(-3) 

4.0 

3.319 

3.393(-l) 

3.595(-l) 

-3.653(-l) 

1.103(-2) 

3.8 

3.179 

3.623(-l) 

3.520(-l) 

-3.680(-l) 

1.378(-2) 

3.6 

3.042 

3.874(-l) 

3.447(-1) 

-3.707(-l) 

1.731(-2) 

3.4 

2.910 

4.145(-1) 

3.376(-l) 

-3.736(-l) 

2.183(-2) 

3.2 

2.786 

4.434(-l) 

3.309(-l) 

-3.765(-l) 

2.758(-2) 

3.0 

2.670 

4.737(-1) 

3.250(-l) 

-3.792(-l) 

3.481 (-2) 

2.9 

2.617 

4.890(-l) 

3.223(-l) 

-3.805(-l) 

3.902(-2) 

2.8 

2.567 

5.042(-l) 

3.200(-l) 

-3.817(-1) 

4.361(-2) 

2.7 

2.521 

5.192(-1) 

3.181(-1) 

-3.827(-1) 

4.860(-2) 

2.6 

2.480 

5.334(-l) 

3.166(-1) 

-3.835(-l) 

5.384(-2) 

2.5 

2.445 

5.467(-l) 

3.156(-1) 

-3.840(-l) 

5.927(-2) 

2.4 

2.415 

5.587(-l) 

3.153(-1) 

-3.842(-l) 

6.475(-2) 

2.3 

2.391 

5.691(-1) 

3.156(-1) 

-3.840(-l) 

7.008(-2) 

2.2 

2.375 

5.774(-l) 

3.165(-1) 

-3.835(-l) 

7.491(-2) 

2.1 

2.366 

5.830(-l) 

3.180(-1) 

-3.827(-1) 

7.874(-2) 

2.0 

2.364 

5.853(-l) 

3.195(-1) 

-3.818(-1) 

8.086(-2) 


Table 5. Equilibrium sequence of irrotational Roche—Riemann 
type binary systems with AT|_/M_=0.5 


Figure 2. Same as Figure |lj but for quantities of a sequence of 
Roche type binary systems with M+/M_ = 0.1. 


alescence. Our numerically exact results for an irrotational 
Darwin-Riemann sequence with equal mass are tabulated in 
Table fij. In Figures g](a)-(c) we show our results and those 
of LRS2. Two results agree with each other very well even 
near the turning point of the sequence to within 0.5%. Two 
fluid components of our present computation come to con¬ 
tact each other at the point of the smallest separation in 
Figure |E] and Table [l]. 

The turning point of J or E corresponds to the dy¬ 
namical instability point of the sequence just as the Roche- 
Riemann binary systems. The dynamical instability point 
appears before the two stars come to contact as the sepa¬ 
ration decreases. This suggests that the binary configura¬ 
tion becomes dynamically unstable when the separation de¬ 
creases due to gravitational radiation emission before these 
two fluid components contact each other. This behavior of 
the sequence is the same as that obtained in LRS2. The sep¬ 
aration of the sequence where J or E attains its minimum is 
larger by 6% and the corresponding frequency is 10% smaller 


d 

da 

n 

J 

E 

dE/dt 

5.0 

3.619 

2.977(-l) 

3.334(-l) 

-2.046C-1) 

5.607(-3) 

4.8 

3.495 

3.139(-1) 

3.277(-l) 

-2.063(-l) 

6.700(-3) 

4.6 

3.372 

3.313(-1) 

3.220(-l) 

-2.081(-1) 

8.037(-3) 

4.4 

3.252 

3.500(-l) 

3.164(-1) 

-2.101 (-1) 

9.674(-3) 

4.2 

3.135 

3.700(-l) 

3.109(-1) 

-2.121 (-1) 

1.168(-2) 

4.0 

3.022 

3.913(-1) 

3.055(-l) 

-2.141(-1) 

1.413(-2) 

3.8 

2.914 

4.137(-1) 

3.004(-l) 

-2.162(-1) 

1.710(-2) 

3.6 

2.812 

4.371(-1) 

2.956(-l) 

-2.183(-1) 

2.067(-2) 

3.4 

2.719 

4.608(-l) 

2.913(-1) 

-2.202(-l) 

2.488(-2) 

3.2 

2.634 

4.843(-l) 

2.876(-l) 

-2.219(-1) 

2.971(-2) 

3.0 

2.563 

5.064(-l) 

2.849(-l) 

-2.233(-l) 

3.499(-2) 

2.9 

2.533 

5.165(-1) 

2.840(-l) 

-2.238(-l) 

3.769(-2) 

2.8 

2.511 

5.240(-1) 

2.835(-l) 

-2.240(-l) 

3.984(-2) 

2.7 

2.492 

5.312(-1) 

2.833(-l) 

-2.242(-l) 

4.208(-2) 

2.6 

2.473 

5.385(-l) 

2.833(-l) 

-2.241 (-1) 

4.452(-2) 

2.5 

2.458 

5.457(-l) 

2.839(-l) 

-2.238(-l) 

4.725(-2) 

2.4 

2.452 

5.489(-l) 

2.848(-l) 

-2.234(-l) 

4.869(-2) 

2.3 

2.453 

5.499(-l) 

2.856(-l) 

-2.229(-l) 

4.934(-2) 


than those of LRS’s result. These differences from those of 
LRS are larger than those for the point mass-fluid star bi¬ 
nary systems. It is because, for the fluid star-fluid star binary 
system, both components are subject to tidal effect from the 
other companion star. 

In Table |s], we tabulate a sequence of irrotational fluid- 
fluid binary configurations with M+/M- = 0.5. Physical 
quantities of this sequence are also plotted in Figure 
By considering that the mass of neutron stars ranges from 
the minimum mass of ~ 0.7 Mq to the maximum mass of 
~ 1.4Mq, it is reasonable to take the minimum of the mass 
ratio to M+/M- = 0.5. For any mass ratio, our results co¬ 
incide qualitatively with those of LRS again, except for the 
sequences of the mass ratios around M+/M- =0.1. Since 
the primary star deforms considerably for models on such 
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Table 6. Equilibrium sequence of irrotational Roche-Riemann 
type binary systems with Mf|_/M_=0.1 


d 

d-G 


J 

E 

dE/dt 

10.0 

4.614 

2.066(-l) 

1.398(-1) 

-3.226(-2) 

2.280(-4) 

9.6 

4.444 

2.186(-1) 

1.373(-1) 

-3.282(-2) 

2.752(-4) 

9.2 

4.275 

2.316(-1) 

1.346(-1) 

-3.341(-2) 

3.341(-4) 

8.8 

4.108 

2.460(-l) 

1.320(-1) 

-3.405(-2) 

4.082(-4) 

8.4 

3.943 

2.616(-1) 

1.293(-1) 

-3.473(-2) 

5.015(-4) 

8.0 

3.780 

2.787(-l) 

1.266(-1) 

-3.545(-2) 

6.200(-4) 

7.6 

3.620 

2.975(-l) 

1.239(-1) 

-3.623(-2) 

7.707(-4) 

7.2 

3.463 

3.179(-1) 

1.213(-1) 

-3.706(-2) 

9.631(-4) 

6.8 

3.311 

3.402(-l) 

1.186(-1) 

-3.793(-2) 

1.208(-3) 

6.4 

3.165 

3.642(-l) 

1.160(-1) 

-3.885(-2) 

1.520(-3) 

6.0 

3.026 

3.898(-l) 

1.135(-1) 

-3.980(-2) 

1.912(-3) 

5.6 

2.896 

4.167(-1) 

l.lll(-l) 

-4.076(-2) 

2.395(-3) 

5.2 

2.779 

4.439(-l) 

1.090(-1) 

-4.168(-2) 

2.973(-3) 

4.8 

2.677 

4.703(-l) 

1.072(-1) 

-4.250(-2) 

3.627(-3) 

4.4 

2.596 

4.937(-l) 

1.059(-1) 

-4.314(-2) 

4.304(-3) 

4.3 

2.579 

4.987(-l) 

1.057(-1) 

-4.326(-2) 

4.465(-3) 

4.2 

2.565 

5.033(-l) 

1.055(-1) 

-4.336(-2) 

4.616(-3) 

4.1 

2.552 

5.074(-l) 

1.053(-1) 

-4.343(-2) 

4.760(-3) 

4.0 

2.542 

5.110(-1) 

1.052(-1) 

-4.347(-2) 

4.888(-3) 

3.9 

2.533 

5.140(-1) 

1.052(-1) 

-4.349(-2) 

5.001(-3) 

3.8 

2.528 

5.162(-1) 

1.052(-1) 

-4.349(-2) 

5.092(-3) 

3.7 

2.524 

5.177(-1) 

1.053(-1) 

-4.345(-2) 

5.158(-3) 

3.6 

2.523 

5.184(-1) 

1.054(-1) 

-4.339(-2) 

5.199(-3) 

Table 7. Equilibrium sequence of irrotational Darwin-Riemann 

type 

binary systems with 

equal mass 



d 

do 

n 

J 

E 

dE/dt 

5.0 

4.060 

2.506(-l) 

3.973(-l) 

-6.589(-l) 

4.002(-3) 

4.8 

3.910 

2.653(-l) 

3.900(-l) 

-6.608(-l) 

4.848(-3) 

4.6 

3.761 

2.814(-1) 

3.826(-l) 

-6.628(-l) 

5.911(-3) 

4.4 

3.613 

2.989(-l) 

3.753(-l) 

-6.649(-l) 

7.254(-3) 

4.2 

3.468 

3.181(-1) 

3.679(-l) 

-6.672(-l) 

8.965 (-3) 

4.0 

3.326 

3.391(-1) 

3.606(-l) 

-6.696(-l) 

l.U5(-2) 

3.8 

3.187 

3.620(-l) 

3.535(-l) 

-6.721(-1) 

1.397(-2) 

3.6 

3.053 

3.869(-l) 

3.466(-l) 

-6.747(-l) 

1.761(-2) 

3.4 

2.925 

4.138(-1) 

3.402(-l) 

-6.773(-l) 

2.230(-2) 

3.2 

2.805 

4.422(-l) 

3.345(-l) 

-6.798(-l) 

2.833(-2) 

3.1 

2.749 

4.568(-l) 

3.320(-l) 

-6.809(-l) 

3.194(-2) 

3.0 

2.696 

4.714(-1) 

3.298(-l) 

-6.819(-1) 

3.597(-2) 

2.9 

2.648 

4.861(-1) 

3.281(-1) 

-6.828(-l) 

4.045(-2) 

2.8 

2.604 

5.002(-l) 

3.270(-l) 

-6.834(-l) 

4.533(-2) 

2.7 

2.565 

5.138(-1) 

3.264(-l) 

-6.837(-1) 

5.058(-2) 

2.6 

2.533 

5.264(-l) 

3.265(-l) 

-6.836(-l) 

5.616(-2) 

2.5 

2.505 

5.385(-l) 

3.275(-l) 

-6.831(-1) 

6.228(-2) 

2.4 

2.489 

5.468(-l) 

3.291(-1) 

-6.823(-l) 

6.716(-2) 

2.3 

2.478 

5.538(-l) 

3.315(-1) 

-6.810(-1) 

7.202(-2) 

2.2 

2.475 

5.581(-1) 

3.341(-1) 

-6.795(-l) 

7.573(-2) 

2.1 

2.478 

5.596(-l) 

3.365(-l) 

-6.782(-l) 

7.775(-2) 

2.0 

2.482 

5.588(-l) 

3.371(-1) 

-6.778(-l) 

7.755(-2) 


sequences, results of our computations become less accurate 
as far as we employ the grid numbers mentioned before. 


4 DISCUSSION AND CONCLUSION 

Since we have succeeded in developing a highly accurate 
numerical code, we can investigate fully deformed configura¬ 
tions for various types of incompressible fluid binary systems 
as far as configurations with the ^-independent planar flows 
are concerned. Among them we have concentrated on se- 


Figure 3. Same as Figure jlj but for quantities of a sequence 
of irrotational Roche-Riemann type binary systems with equal 
mass. 


quences of the Roche, the irrotational Roche-Riemann and 
the irrotational Darwin-Riemann type configurations in this 
paper because these sequences are physically important for 
the evolution of BH-NS or NS-NS binary systems due to 
gravitational radiation emission. Concerning general binary 
systems for which the vorticity seen from the inertial frame 
does not vanish and the internal motion exists in the ro¬ 
tating frame, such as Roche-Riemann or Darwin-Riemann 
binaries, we have computed several sequences. However, if 
the initial spin of the fluid star is rather small compared 
with the final orbital angular velocity, the spin does not af¬ 
fect the final results appreciably. Therefore we did not show 
our results for such sequences in this paper. 

Differences between our results and those of LRS be¬ 
come relatively large when the two component stars are 
coming closer and the mass ratio of two components differs 
from unity. It is because stars deform considerably around 
the Roche lobe overflowing stage which cannot be expressed 
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Figure 4. Same as Figure H but for M+/M- = 0.1. 


by the ellipsoidal approximation used by LRS. However the 
qualitative character of the equilibrium sequences are the 
same. In other words, the dynamical and/or secular insta¬ 
bility point appears before the models reach the Roche lobe 
overflowing stage, as shown in LRS1 and LRS3. Our results 
show that generally the separation at the onset of the in¬ 
stability is larger and correspondingly the orbital frequency 
at that point is smaller compared to the results of the ellip¬ 
soidal approximation by the amount of 2 ~ 10%. 

In order to clarify the final fate of binary systems, 
it would be better to define four critical distances as fol¬ 
lows (see e.g. Lai, Rasio & Shapiro 1993b). (i) tr\ the 
Roche (Darwin)-Riemann limit below which no equilibrium 
states exist, i.e. the minimum separation of binary stars; 
(ii) rc- the contact limit where two stars contact each other 
(note that this limit exists only for a fluid-fluid star type bi¬ 
nary system with equal mass and the equilibrium sequence 
may continue to a Dumbbell type sequence); (iii) r m : the 
hydrodynamical stability limit below which stars fall down 
onto each other dynamically due to the tidal potential; and 


Figure 5. Same as Figure jl| but for quantities of a sequence 
of irrotational Darwin-Riemann type binary systems with equal 
mass. 


(iv) vgr- the general relativistic stability limit where no sta¬ 
ble orbits exist due to the effect of strong gravity. Although 
the last critical distance vgr does not appear in our present 
treatment, the binary stars fall down onto each other when 
two stars come within tgr ~ 6GMtot/c 2 . 

Concerning a typical BH-NS binary system whose 
masses are around ~ 10Mg (BH) and ~ lJU© (NS), we 
should consider the general relativistic effect of BH at the 
final stage of evolution. Taniguchi & Nakamura ( 199f| ) have 
computed the incompressible ellipsoidal model of an irrota¬ 
tional neutron star around a black hole which is represented 
by using the improved pseudo-Newtonian potential. In their 
paper they studied models for which relativistic effect domi¬ 
nates, i.e. the models with tr < r m < vgr as typical BH NS 
systems. It implies that such binary systems become dynam¬ 
ically unstable due to the general relativistic effect but not 
due to the tidal effect and that the dynamical instability 
occurs before the star deforms significantly. 
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Table 8. Equilibrium sequence of irrotational Darwin—Riemann 
type binary systems with My /M— =0.5 


d 

d-G 

D 

J 

E 

dE/dt 

5.0 

4.420 

2.205(-l) 

3.683(-l) 

-6.875(-l) 

2.056(-3) 

4.8 

4.250 

2.339(-l) 

3.612(-1) 

-6.892(-l) 

2.506(-3) 

4.6 

4.081 

2.486(-l) 

3.540(-l) 

-6.909(-l) 

3.078(-3) 

4.4 

3.913 

2.650(-l) 

3.467(-l) 

-6.928(-l) 

3.812(-3) 

4.2 

3.746 

2.830(-l) 

3.394(-l) 

-6.948(-l) 

4.761(-3) 

4.0 

3.580 

3.031(-1) 

3.320(-l) 

-6.970(-l) 

6.002(-3) 

3.8 

3.417 

3.254(-l) 

3.246(-l) 

-6.993(-l) 

7.640(-3) 

3.6 

3.256 

3.502(-l) 

3.172(-1) 

-7.018(-1) 

9.825(-3) 

3.4 

3.099 

3.778(-l) 

3.100(-1) 

-7.044(-l) 

1.277(-2) 

3.2 

2.948 

4.083(-l) 

3.031(-1) 

-7.071(-1) 

1.677(-2) 

3.0 

2.804 

4.417(-1) 

2.969(-l) 

-7.098(-l) 

2.222(-2) 

2.8 

2.674 

4.773(-l) 

2.919(-1) 

-7.121(-1) 

2.963(-2) 

2.75 

2.644 

4.864(-l) 

2.910(-1) 

-7.126(-1) 

3.186(-2) 

2.7 

2.615 

4.953(-l) 

2.902(-l) 

-7.130(-1) 

3.422(-2) 

2.65 

2.589 

5.042(-l) 

2.897(-l) 

-7.133(-1) 

3.676(-2) 

2.6 

2.565 

5.128(-1) 

2.894(-l) 

-7.134(-1) 

3.938(-2) 

2.55 

2.543 

5.211(-1) 

2.894(-l) 

-7.134(-1) 

4.218(-2) 

2.5 

2.525 

5.285(-l) 

2.898(-l) 

-7.132(-1) 

4.495(-2) 

2.45 

2.511 

5.351(-1) 

2.907(-l) 

-7.128(-1) 

4.769(-2) 

2.4 

2.503 

5.404(-l) 

2.924(-l) 

-7.119(-1) 

5.036(-2) 


Since the point of dynamical instability always appears 
prior to the Roche (Darwin)-Riemann overflowing state as 
shown by LRS as far as the ellipsoidal approximation is con¬ 
cerned, a possibility of occurrence of the Roche (Darwin)- 
Riemann overflow has not been considered seriously. There¬ 
fore their conclusion is simple, i.e. that binary systems nec¬ 
essarily result in dynamical coalescence, although general 
relativity, tidal effect and compressibility may change the 
results. On the other hand, for realistic BH-NS binary sys¬ 
tems, there may exist configurations for which the Roche- 
Riemann overflow will occur prior to dynamically falling as 
an alternative final fate, i.e. models with r m < tr. In order 
to show whether such a situation can be realized, we should 
take the compressibility of the fluid component into consider¬ 
ation because we have shown that there always exist dynam¬ 
ical instability points on the incompressible BH-NS binary 
sequences. We can expect that the effect of compressibility 
leads the fluid star to reach the Roche-Riemann overflow 
state before getting to a dynamically unstable state. It is 
opposite to the effect of strong gravity of BH. This effect 
of compressibility will be discussed in the subsequent paper 
(Uryu & Eriguchi 1998a). 

If the equilibrium sequence of a binary system cannot 
reach a point where dynamical instability sets in, we need 
to consider the problem whether this Roche-Riemann over¬ 
flow proceeds violently or not. It will certainly depend on 
the compressibility and the gravitational field of the black 
hole. Roughly speaking, for the fluids with stiff equations of 
state, if the star loses mass, the radius of the star is likely 
to shrink. Thus the mass overflow will stop. However, since 
there is gravitational wave emission, the system will evolve 
to the Roche-Riemann overflow stage again. In this case the 
mass overflow will continue on the time scale of gravitational 
wave emission. On the other hand, if the equation of state 
is soft, the radius of a mass losing star will expand. It im¬ 
plies that the mass overflow will continue violently. Such a 
mass transfer between two stars has been discussed by many 
authors (see Kochanek 1992; Bildsten & Cutler 1992 and 
references therein). When we consider the black hole, this 


Figure 6. Same as Figure H but for M_|_/Af_ = 0.5. 


situation may be analogous to the runaway mass overflow 
for self-gravitating axisymmetric toroidal configurations (see 
e.g. Nishida et al. 1996; Masuda & Eriguchi 1997). 

Although we have only considered incompressible 
fluids in this paper, the compressibility would make 
Roche (Darwin)-Riemann overflow more probable. Conse¬ 
quently the final fate of real binary systems with non-equal 
masses could be determined from the competition between 
the occurrence of the state of the marginally stable orbit 
due to general relativity and the occurrence of the state of 
Roche (Darwin)-Riemann overflow. For compressible gases, 
we have recently succeeded in developing a new computa¬ 
tional metho d for irrotational gaseou s binary systems with 
equal mass (Uryu & Eriguchi 1998a). We will be able to 
investigate more realistic situations in the near future as 
long as Newtonian gravity is concerned, and the results pre¬ 
sented here will be extensively used to calibrate the results 
computed by the new computational code. Concerning the 
distance rGR, we will be able to implement the general rel¬ 
ativistic effect in our treatment of the Roche-Riemann type 
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binary systems by using the pseudo-Newtonian potential as 
was done by Taniguchi & Nakamura (1996). After such com¬ 
putations, the nature of the final fate will be finally clarified. 
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APPENDIX A: BASIC EQUATIONS FOR 
ACTUAL COMPUTATIONS 

In this section, we describe the basic equations used in our 
actual computations. To clarify the discussion hereafter, we 
call the stellar surface R(9, p) and Fourier coefficients a m as 
the physical variables or the variables and the other quan¬ 
tities, £, 12, d and C as the physical parameters or the con¬ 
stants. They are the unknowns which appear in our formu¬ 
lation as will be seen below. 

Because of the incompressibility and the choice of the 
coordinate mentioned in section 2, we can integrate the grav¬ 
itational potential of the fluid and express it in terms of the 


stellar surface R(9, p). It means that the value of the grav¬ 
itational potential of the incompressible fluid is determined 
only from the shape of the stellar surface R{9,p). Conse¬ 
quently we can reduce the physical variables in the Bernoulli 
equation (jiij) to the stellar surface R{9 , p), the Fourier coef¬ 
ficients o m and other constants. Therefore we will describe 
each term of this Bernoulli equation on the surface. 

The first term of the LHS of equation (0) on the sur¬ 
face is calculated after substituting equation ( jll| ) into ([!]) as 
follows: 

u x + u y ) = ^C 2 |^ [R(0,V>) sin<9] 2 


+ ^ ' m am [R(9,p) sin#] m cosmip 

ri=l 
oo 

'y ' ma m [R(9,p) sin0] m-1 sin m p 


+ 


'y ' m am [R(9,p) sin#]™ 1 cos mp 


(Al) 


The second term of the LHS of equation (|2j) on the 
surface is directly expressed from equation (ju|) as 

(C + = -(C + 2Q.)A-\R{9,p) sin 9] 2 


OO N 

+ y ' a m [R(9,p) sin#] m cos mip > . (A2) 

™-n ' 


For the fluid star-fluid star binary system, we need to add 
subscripts ± in the two equations shown above. The third 
term of the LHS of equation (j2|) on the surface is calculated 
as follows for the fluid star-fluid star binary: 

-\tf{X 2 ± + Yl) = {[R±(9±,p±) sin 9±} 2 

± 2d± R±(9±, p±) sin0± cos<p± + d±} . (A3) 

The fourth term of the LHS vanishes identically on the stel¬ 
lar surface. 

The fifth term of the LHS of equation (^|) which is the 
gravitational potential term is divided into two parts: 1) the 
contribution from its own component (f> s and 2) that from 
the other (j> £ , i.e. the relation (f> = <j> 3 + (j> e holds. These grav¬ 
itational potentials are expressed in the integral form of the 
Poisson equation (see equation (f^)). Since the Green’s func¬ 
tion can be expanded in terms of the Legendre functions, 
we can analytically integrate the potentials along the ra¬ 
dial direction (see e.g. Eriguchi, Hachisu & Sugimoto 1982; 
Eriguchi & Hachisu 1983; Eriguchi & Hachisu 1985). After 
the integration, the gravitational potentials on the stellar 
surface can be expressed as follows. The self-gravity of its 
own component becomes 

/*7T p2tT 

4> s (6, ip) = — Gp / sin 0'd9 r / dtp' 

Jo Jo 


X Y, fn [R(9 , p), R(0’, p')] Pn (cos (3), (A4) 

n =0 

where f n [R(9,p), R(9',p')\ is defined as 
f n [R(9,p),R(9',p')] = 
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' 1 R { e', v >') n+3 

n + 3 R(9,ip) n+1 
for R(0',<p) <R{0,y), 

(—— + —!—) R(0,ip) 2 

In + 3 n — 2/ y 

1 R(Q,<p) n 

n- 2 R(6',ip') n - 2 
for R(9, tp) < R(9' , 7 /) and n ^ 2, 


(A5) 


P ( 0,¥>) 2 


1 _RM_ 

5 -R( 6 »',^) 


, for R(9,ip) < R(9',ip') and n = 2. 

Here P n is the Legendre function and cos /3 is defined as 
cos/3 = cos 9 cos 0* + sin 9 sin 9' cos(tp — ip'). (A 6 ) 


In the above expression we omit subscripts ± which should 
be used to distinguish the primary and the secondary stars. 

The external gravitational potential from the other 
component can be expressed as follows: 


<f>e±(9±,<p±) = —Gp sin 9 T d9^ 




E 




n + 3 


D 


n-\-1 


■ P„(cos 7 ±), 


(A7) 


where D± and cos 7 ± are defined as 
D± = {( d++d _) 2 + R±(9±,p±) 2 

±2 (d+ + d-)R±(9±, ip±) sin#± cosip±} 1 ^ 2 , (A 8 ) 


and 

COS7 = (±(d+ + d-)R^(9^, p T ) sin# T cos<p T 
+ R±{9±, p±)R^(9^, tpY) cos /3} /R^- (^=fj <p^)D±, (A9) 

respectively. Because of the symmetry of the x-z and the x-y 
planes, we can impose the following condition for the shape 
of the star: 


R(9, p) = R(n — 9, <p) = R(9, 2n — p). (A10) 


Consequently we can reduce the integral region to a quarter 
of the whole space, 



(All) 


Accordingly, the Legendre polynomials in the gravitational 
potentials should be rearranged. For P n (cos/3) in <j> a , we 
should replace it as follows: 


Pn (cos 0) -> Pn (COS (3a ) + Pn (COS /3b) 

+ P„(cos/3c) + P„(cos/3d), (A 12 ) 

where 

cos (3 a = cos 9 cos 9' + sin 9 sin 9' cos(ip — p), 

cos fib = — cos 9 cos 9' + sin 9 sin 9' cos(<p — p'), 

COS (5 C = <X>s9cOs9' +SV0.9sm.9'COs{p-\-p), 

cos Pd = — cos 9 cos 9' + sin 0 sin 0' cos(t> + tp'). (A13) 


The quantity P„( cos 7 ) in the external potential is rear¬ 
ranged in the same way. 

As we mentioned in section 2, the Bernoulli equation is 
discretized on the surface grid points. The discrete approxi¬ 
mation to the Bernoulli equation gives one equation for each 
point on the surface grid and can therefore be used 

to determine the surface R(9i g ,pi v ). 

We have another physical variable a m in our formula¬ 
tion. As mentioned in section 2, they are determined from 
the boundary condition of the stream function 'L. By using 
equation (|ll|), equation (^) can be expressed as follows: 

OO 

-R(n/2,p) 2 + amR(n/2, p)" 1 cos mp = 0. (A14) 

m=0 

This equation is also discretized on the grid points of pi v . 
Since we include the coefficients up to 0 < m < N v / 2, if we 
adopt the grid points with even integers of i v on the equa¬ 
torial plane, the number of the equations and the number of 
the coefficients coincide. As an alternative way, it is possible 
to use all the grid points on the equatorial surface of the 
star and employ the least square method to determine o m . 


APPENDIX B: CONDITIONS FOR 
PARAMETERS 


In this section we describe the scheme to determine the phys¬ 
ical parameters or the constants. The choice of the param¬ 
eters is essential to construct a successful numerical code 
which can be used to get converged equilibrium solutions. 
In our actual computations, we use non-dimensional vari¬ 
ables as follows: 



Ci = 

c = 

C rn, “ 


n 

YGP 

c 

GpR 2 0 ’ 


Pn 


Ui 

~ YGpRo’ 


T 

YG~pRl' 




s YG-p 

V 

p = —, 

Pc 

3 = Pc 
P Gp 2 R 2 0 ’ 

y = E (bi) 


where Ro is the half of the geometrical width of the primary 
star along the x axis, i.e. Ro = (P+ ut — R+)/2 and p c is the 
pressure at the coordinate center of the each component. 

Using these non-dimensional variables, we can rewrite 
the Bernoulli equation (j2[) as follows: 

^(fA + Wy) + (C + 2fi)'F 

- if) 2 (X 2 + Y 2 ) + ftp + 4> = C. (B2) 


We note that as we have shown in the previous section 
we only need to handle the equations on the stellar sur¬ 
face where the pressure vanishes. Therefore the parameter 
/3 does not appear in the equation but can be determined 
afterwards. 

For the fluid star-fluid star binary system, there are 
seven physical parameters, i.e. 


C±> d± ,C±. 


(B3) 


Therefore we should impose seven conditions so as to de- 
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termine these parameters. As for the quantities £±, we can 
specify some simple relations of £ and Cl as 

F±(C±,0) = 0. (B4) 

For the circulation preserving sequences, this relation should 
be 

C± + Cl = Co±, (B5) 

where £o± are given constants. For the synchronously rotat¬ 
ing sequences, we can choose the relation 

C± = 0. (B 6 ) 

As for the distances from the rotational axis to the origins of 
the spherical coordinate systems of each component d± , two 
conditions are required. The total separation d+ + d- and 
the mass ratio M+/M- are specified for this purpose. Since 
the origins of the spherical coordinates of two component 
stars are chosen to be located at the geometrical centers of 
the two stars on the x axis, we can impose the following 
three conditions: 


J 4 ( 7 t/ 2 , 0 ) = 1, R+( 7 r/ 2 , 7 r) = l, 

R-(n/2,0) = R-(tv/2,tv). (B7) 


It should be noted that the physically meaningful conditions 
to specify the parameters are four, in other words, they are: 
1 ) two conditions for £±, 2 ) one condition for the total sepa¬ 
ration and 3) one condition for the mass ratio of the binary 
system. 

For the point source-fluid star binary system, the phys¬ 
ical constants related to the secondary point source are dif¬ 
ferent from those of the fluid-fluid binary system. In this 
case parameters and C_ appear no more but we need 
to include one parameter which characterizes the gravity of 
the secondary star. In this paper we choose the mass of the 
secondary as a parameter: 


(ps+ 


GM _ 

“dT" 


(B 8 ) 


Therefore the total number of parameters becomes six for 
this case: 


C+, n, d± ,C+, M- 


(B9) 


Consequently we should specify six relations for them, five 
of which are the same as the previous case: we need not 
use the relation for and the relation for the radius of 
the secondary star. The sixth relation is obtained from the 
requirement that the motion of the point mass should be 
Keplerian subjected to the gravitational force of the pri¬ 
mary star. It can be written by using the non-dimensional 
variables as follows: 


d-Cl 2 


OO 


E 

n =0 


(d+ + d_)»+ 2 ’ 


(BIO) 


where 


In = 


n + 1 
n + 3 


sin 9+dd 4 


dip+ 


x R + (6 +, ifi+) P„ (cos 7 ). 


(Bll) 


Therefore, the physically meaningful conditions for the pa¬ 
rameters are three: 1 ) the relation for £+, 2 ) one condition 


for the total separation and 3) one condition for the mass 
ratio of the binary system. 

These equations for the parameters are solved simul¬ 
taneously with the physical variables mentioned before by 
using the Newton-Raphson scheme. This method is so pow¬ 
erful that we can obtain converged solutions after only 5 or 
6 iteration cycles. 
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